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We study the non-degenerate optical parametric oscillator in a planar interferometer near thresh¬ 
old, where critical phenomena are expected. These phenomena are associated with non-equilibrium 
quantum dynamics that are known to lead to quadrature entanglement and squeezing in the oscil¬ 
lator field modes. We obtain a universal form for the equation describing this system, which allows 
a comparison with other phase transitions. We find that the unsqueezed quadratures of this system 
correspond to a two-dimensional XY type model with a tricritical Lifshitz point. This leaves open 
the possibility of a controlled experimental investigation into this unusual class of statistical models. 

We evaluate the correlations of the unsqueezed quadrature using both an exact numerical simulation 
and a Gaussian approximation, and obtain an accurate numerical calculation of the non-Gaussian 
correlations. 


I. INTRODUCTION 

Non-equilibrium pattern-formation occurs in many 
physical systems, giviim rise to the emergence of order 
on macroscopic scales [l| . The theory of hydrodynamics 
is a paradigm for understanding these phenomena, and 
is applicable to many branches of physics, chemistry, bi¬ 
ology, astrophysics, and other sciences 0. This is often 
applied to many body systems subject to nonlinear cou¬ 
pling in a dissipative environment with external fluxes. 
In physics, one of the most studied hydrodynamic effects 
is the theory of fluid flows near the Rayleigh-Benard in¬ 
stability instability [l|. The simplest nontrivial model 
for fluid convection that displays pattern formation due 
to convective instabilities was derived by Swift and Ho- 
henberg , where nonlinear coupling of fluctuations was 
included to demonstrate the failure of mean field theory 
for critical exponents. 

Here, we treat pattern formation and universality in a 
paradigmatic non-equilibrium quantum system: the non¬ 
degenerate parametric oscillator (OPO). Theories of a 
similar nature have been applied to non-equilibrium spa¬ 
tially extended structures in lasers and other related sys¬ 
tems, with an emphasis on universal behavior of phase- 
transitions, pattern formation and self-organization Q. 
Yet down-conversion in a non-degenerate parametric sys¬ 
tem can display new possibilities not found in these sim¬ 
pler cases. In particular, one can have entanglement and 
EPR paradoxes ii. 

In the present paper, we investigate a new type of 
critical point phase transition in this non-equilibrium 
quantum system in order to understand the universality 
class. In a type II OPO, there are two down-converted 
fields with orthogonal polarization. Hence, the order pa¬ 
rameter is a complex or vector field in two dimensions. 
The two components of this vector field are associated 
with the polarization degrees of freedom of the down- 
converted radiation field. This is a quantum system 
driven to a phase transition far from thermal equilibrium. 


It is also known to display strong quantum entanglement 
and EPR correlations Q in the case where there are two 
correlated output modes. We wish to understand the 
behavior of this phase-transition using a first-principles 
analysis 0,0 of the relevant master equation. 

Experimentally, the OPO is now a mature technol¬ 
ogy with both commercial and fundament al ap plications. 
Following initial theoretical predictions the quan¬ 

tum limited type I OPO was investigated experimen¬ 
tally by Wu et al [a, demonstrating quantum squeez¬ 
ing. Later the type H case, in a triply resonant cav¬ 
ity, was used to experimentally demonstrate continuous 
variable EPR correlations 0, also originally predicted 
theoretically Bi- These initial experimental investiga¬ 
tions were in few mode devices. Both below and above 
threshold experiments have been carried out, as close as 
±1% of the critical point [iM3> confirming predictions 
of thresholds and conversion efficiency. Operation at the 
critical point results in low-frequency critical fluctuations 
and non-Gaussian behavior [iSj . 

Spatially extended pattern formation in a degenerate 
or type I OPO has also been analyzed previously [il. It 
is related to the Lifshitz phase transition This is a 
model used to describe the phase transition to a modu¬ 
lated magnetic phase [2l|, . In this simplest case one 

has a two dimensional, planar system with a scalar order 
parameter , which has known universality properties. 
A Swift Hohenberg equation was derived for spatially ex¬ 
tended nondegenerate type I OPOs with flat end mirrors 
[IJ, but ignoring fluctnations. OPO experiments have 
been recently extended to these types of multi-mode de¬ 
vices [ 2 ^, [2g, with type H as well as type I paramet¬ 
ric downconversion. The experimental situation is that 
while multimode experiments have mostly used confocal 
mirrors, the first type H multimode planar mirror exper¬ 
iments have been carried out [ 2 ^. 

Studies of the Swift Hohenberg equation in the neigh¬ 
borhood of the critical point, the Lifshitz point, as well 
as pattern formation for lasers have been discussed 
The Lifshitz point is similar to the tricritical point 281. 
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Tricritical points occur in different physical systems. 
They correspond to a point in a phase-diagram where 
two lines of ordinary critical points meet and terminate 
(^ . For critical points and tricritical points there are 
two critical dimensions. The upper critical dimension 
refers to the one above which the critical exponents have 
classical values (2^. The lower critical dimension refers 
to the smallest dimension for which there is a true phase- 
transition, due to increasing fluctuations as the dimen¬ 
sion is decreased. 

Systems like these can entangle large numbers of 
modes, with quantum entanglement increasing as thresh¬ 
old is approached. Currently, experiments near the crit¬ 
ical point are sensitive to classical fluctuations [l^ and 
heating effects [2^. With improved stabilization meth¬ 
ods, we expect that these technical problems can be over¬ 
come. However, noise due to quantum fluctuation effects 
will remain. 

Theoretically, the usual approach for two-dimensional 
non-equilibrium problems is the Landau-Ginzburg (LG) 
equation (^ . which was first used for understanding su¬ 
perconductivity near threshold. Normally, the LG equa¬ 
tion is derived using approximations like adiabatic elimi¬ 
nation, and is taken as an effective equation for the phys¬ 
ical system. There is now an increased interest in ex¬ 
tended spatial and multi-mode structures HI,!!! in the 
quantum optical parametric oscillator (OPO) 'um. 
due to availability of multiple transverse mode cavi¬ 
ties (33| and quantum imaging control [^ . It is im¬ 
portant in these cases to understand how quantum noise 
enters the dynamical equations (^ . 

We show that planar type-II down-conversion creates 
a non-equilibrium system with a similar general type of 
symmetry to the Berezinskii and Kosterlitz-Thouless 
(BKT) model [s^, but with an isotropic Lifshitz point, 
first studied in magnetic systems by Hornreich et al . 
Physically, this is not a classical fluid or magnetic sys¬ 
tem, but rather is a quantum system, driven into a non¬ 
equilibrium critical point far from thermal equilib¬ 
rium. The non-degenerate planar OPO is fundamentally 
different to the usual BKT model, having a quartic rather 
than quadratic momentum-dependence in the linear re¬ 
sponse function. This places it in a similar category 
to the Swift-Hohenberg model and next-nearest neigh¬ 
bor lattice models. Our model can also display strong 
EPR entanglement and other nonclassical properties in 
addition to the Lifshitz point behavior. The parametric 
model therefore provides a novel path to the investigation 
of unusual classical and quantum noise effects. 

Phase transitions with this general type of symmetry 
are continuous yet break no symmetries. There is also no 
ferromagnetism in this system, as this is prohibited by 
the Mermin-Wagner theorem [^. Berezinskii predicted 
a new type of phase with correlations that decay slowly 
with distance with a power law. The phase transition for 
a ferromagnetic phase is prevented by the appearance of 
vortex and anti-vortex pairs. Many quantum phase tran¬ 
sitions in two dimensions belong to this class. The con¬ 


tinuous version of the XY or Ising models are often used 
to model systems that possess order parameters with a 
symmetry of this type, e.g. superffuid helium, liquid crys¬ 
tals, two dimensional Bose Einstein condensate (BEG), 
and others. While having the same order-parameter and 
space dimension as this case, we show that the type II 
OPO has a different universality class. 

Our approach is to start from the usual master equa¬ 
tion model that describes a type II OPO with transverse 
modes. Here we consider that the OPO is non-degenerate 
in polarization. We then map the coupled Heisenberg 
equations into the positive P-representation of the 
density matrix. This choice is made because the positive 
P-representation allows us to exactly map the density 
matrix evolution into an equivalent stochastic equation. 

We organize this paper in the following way. First 
we describe the Hamiltonian model and derive the equa¬ 
tion of motion for this open system. Next, we use the 
positive-P representation for mapping these equation into 
a Langevin type, which can either be treated numer¬ 
ically or via analytic approximations. We will mostly 
describe the unsqueezed quadratures of the field, which 
have a large similarity with the corresponding magnetic 
system. At this point we can recognize the universality 
class. We give both a Gaussian approximation to the 
correlation functions of the unsqueezed quadratures, and 
high-precision numerical simulations of the non-Gaussian 
corrections. We show that, even though the system is 
far below its upper critical dimension, the non-Gaussian 
character of the fluctuations is relatively small, and the 
intensity fluctuations are nearly factorizable. 

In a following paper we will focus our attention on 
the squeezed field quadratures, to give a spatial map of 
quantum squeezing and Einstein-Podolsky-Rosen (EPR) 
entanglement. 


II. THE MODEL 

The system of interest comprises an optically driven 
planar Fabry-Perot cavity or interferometer with a non¬ 
linear medium that possess a parametric nonlinearity. 
The nonlinear crystal is cut to give a type H phase match¬ 
ing, that couples a pump field to two down-converted 
fields having an orthogonal polarization. The cavity is 
pumped with a spatially extended coherent light with 
frequency wqj with a transverse spatial profile. In the 
simplest case we consider this pump to be a plane wave. 


A. Quantum Hamiltonian 

The outgoing down-converted light, amplified inside 
the cavity, develops structures and patterns due to 
diffraction, nonlinear coupling, and detuning between the 
wavelength of the down-converted field and the cavity 
size. This depends on the modal decomposition of this 
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cavity. The mirrors have parameters that can be con¬ 
trolled by experimentalists. The tunable parameters in¬ 
clude the reflection coefficient for each mode, and the 
cavity detunings for each mode. 

Our model for this system is similar to many earlier 
treatments of driven nonlinear optical cavities p, S, M, 
I 41 I - I 45 I I. It includes a linear coupling between the exter¬ 
nal electromagnetic field modes and the internal cavity 
modes, owing to a partially transmitting mirror. The 
cavity and mirror parameters determine both the cou¬ 
pling to the driving field and the decay rate of the cav¬ 
ity or interferometer. We note that for a low-Q device, 
it is important to use non-orthogonal quasi-modes (d^ . 
Here, we assume the opposite case of a thin, high-Q ex¬ 
tended planar cavity, so that the external modes are sim¬ 
ply plane-wave modes. 

The quantum Hamiltonian in the interaction picture 
has four main terms that can be summarized by the fol¬ 
lowing expression: 

H — Ufree ^int T ^pump ^res : (^■^) 

where there are three fundamental Bose fields, Ai{x,t) 
for i = 0,2. The boson fields are the transverse internal 
modes of the planar cavity, which is assumed to have a 
single longitudinal mode in the direction normal to the 
mirrors. We note that this planar polariton held model 
is also used in the theory of a polaritonic BEC [1^. The 
boson fields have two orthogonal polarizations, * = f,2, 
and obey the usual equal time commutation relation, 


,t) 


(x — x'). 


Here (5^ (x — x') is a Dirac delta function in the two- 
dimensional transverse plane of a; = {x,y). The cav¬ 
ity fields are defined in terms of polaritonic annihilation 
and creation operators as: Ai{x,t) = L, 

where ai{k,t) represents the annihilation operator for a 
free polariton mode with transverse momentum fc, and 
the summation is over a set of M discrete modes with 
periodic boundary conditions, on a large area . 

The free evolution Hamiltonian that accounts for 
diffraction inside the planar cavity is: 


Hfree — 


E*/ 

i=0 


(fx Aj 


Uj „9 


A,,. 


( 2 . 2 ) 


The interaction Hamiltonian representing the coupled 
modes inside a crystal with nonlinearity is given 
by 0: 


Hint — 


I d?x [xioil4 - X*AIA,A 2 \ . (2.3) 


This term may represent a photon of frequency wq con¬ 
verted in two photons of distinct frequencies wi and UJ 2 , 
or photons with orthogonal polarizations, or both. 

From now on, for definiteness, labels f and 2 stand 
for polarizations, and we shall focus on the degenerate 
frequency case with non-degenerate polarization. For di¬ 
mensional reasons, y oc / '/i, where y^^^ is the Bloem- 

bergen nonlinear polarizability coefficient and i is the in¬ 
tracavity longitudinal mirror spacing. 

Following standard input-output theory derivations (^, 
0) , the Hamiltonian term associated with the in¬ 

put laser pumping in a rotating frame at frequency wl 
is: 


Hpump — 


(Px 


r(a;)e2“^‘io - £’(a;)e-2“^*ij 


(2.4) 

While it is possible to choose any shape carrying spatial 
structure for the input pump, here for simplicity we will 
assume a plane wave input. The reservoir Hamiltonian 
is assumed to have the structure: 


Hr 








(2.5) 


i=0 ■ 


Hence, there are local coupling terms to independent ex¬ 
ternal free-field reservoirs for each polarization and each 
spatial mode. The external fields are described by a free 
Hamiltonian 


B. Master equation 

The non-unitary evolution of the system comes from 
the coupling between the cavity modes and the output 
modes. This can be treated as a quantum Markovian 
process that simulates a bath interaction. We carry out 
this calculation in a type of interaction picture so that 
the interaction picture operators evolve according to a 
reference Hamiltonian Hq, given by: 


This Hamiltonian describes a planar cavity with intra¬ 
cavity resonant frequencies uJi and group velocities Vi for 
the three field envelopes. The resonant frequencies have 
the relation loq uii + UJ 2 , where ujq is the fundamen¬ 
tal mode, and uii and uj 2 stand for the down-converted 
light. The case of spherical mirrors can also be analyzed 
in a similar way [33|, but is not treated here. We sup¬ 
press the field space-time arguments to obtain more com¬ 
pact expressions in the integrals. The two-dimensional 
Laplacian is, as usual, jdx^ + d'^ jdy^. A one¬ 

dimensional system can be treated by simply dropping 
one of the dimensions. 


= [ d^xu°A]Aj , ( 2 . 6 ) 

0 

where the reference frequencies ui^ are chosen so that: 

uj° = Wi -I- e Ri wi, 
a ;2 = Wi — e ~ 0J2, 

ujq = 2ujl ~ wq. (2-7) 

We note that while the choice of is determined by the 
pump frequency, the choice of e is arbitrary, as long as 
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the Markovian approximation is still valid. Our choice 
of reference Hamiltonian is not determined by the in¬ 
tracavity frequencies uji, which leads to detuning terms 
occurring in the resulting equations of motion. This al¬ 
lows us to have some freedom of choice in defining the 
detuning parameters. The choice will be made definite 
in the following sections. Following standard techniques 
for deriving the master equations [4l|, |4^ |4^, we can 
write the master equation for the density operator in the 
generalized Lindblad form: 


dp 1 
dt ih . 


H,p 




2=0 


where the dissipative Liouville super-operator, 

[p] = f (fx 2AipAl - pAlAi - AlA^p 


( 2 . 8 ) 


(2.9) 


describes the output coupling of the i-th intracavity mode 
with the external bath. 

We note that, although we focus on the master equa¬ 
tion approach here, it is sometimes useful to write the 
time evolution in a complementary quantum Langevin 
formulation. This formulation is given in Appendix A. 


III. STOCHASTIC EQUATIONS IN THE 
POSITIVE-P REPRESENTATION 

As nonlinear operator equations are not generally sol¬ 
uble, it is more manageable to map the operator equa¬ 
tions into c-number form. In this approach, a master 
equation is transformed into a positive-definite Fokker- 
Planck equation using operator identities that map the 
operator terms in the master equation into differential 
operators 1^ Id^ . To do this we have to use a phase- 

space representation of the master equation. 

Phase-space representations in a classical phase-space 
do not give a positive-definite equation. An example is 
the Wigner representation [d^. While this is exact, it is 
not able to be mapped into Langevin equations, unless 
one either truncates or uses higher order noise [^. One 
can also linearize the Hamiltonian and obtain an approx¬ 
imate Wigner diffusion [s^. Another approach is the 
Husimi Q-function (5ll |. where, in order to obtain a posi¬ 
tive Fokker-Planck equation, an unphysical constraint on 
the phase-space trajectories has to be used (^ . 

Here we wish to have the ability to treat non¬ 
equilibrium structures without restrictions. For this 
purpose, the most useful representation is the positive- 
P representation jd^, which is an extension of the 
Glauber-Sudarshan P-representation into a phase- 
space of double the classical dimensions. Unlike the P- 
representation, which is singular for nonclassical states, 
the positive-P representation is well-defined, positive and 
non-singular for any quantum state. This approach 
allows us to map the density matrix equation into a 
Fokker-Planck equation on a non-classical phase-space. 


In the positive-P representation stochastic averages give 
normally-ordered quantum expectation values. A brief 
description is given in Appendix B. 

The stochastic field partial differential equations are 
given by an extension of our earlier work [23|: 

= —7odlo + £{x) — X* Ai A 2 -I- --2-V^ Ao , 
dt 2 loo 

= -71 Ai -I- x^oAj -I- + V X^oCi) 

= —72A2 + -I- + \/ xAo ^2 -(3.1) 

The three equations that correspond to the hermitian 

conjugate fields, Af, are obtained by conjugating the 
constant terms, and replacing stochastic and noise fields 
so that: Ai —Af and —>■ where and are 

independent Gaussian complex noises. These are equiva¬ 
lent in the mean to the conjugated Heisenberg equations, 
but are independent c-number equations. They are not 
conjugate in every realization. We also note that Ai and 
At are six independent, complex c-number fields. 

The stochastic fields that describe the quantum 
noise are complex and Gaussian, whose non-vanishing 
correlations are: 

{^i{x,t)^ 2 {x',t')) = 5‘^{x - x')S{t - t') 

{^t t)C^ [x', t')) = s^{x - x')5{t - t') . (3.2) 

This means that ^}.{x,t), represent four inde¬ 

pendent, delta-correlated, complex c-number Gaussian 
stochastic fields with zero mean. They are completely 
characterized by the specified correlations. From the 
stochastic equations m, it is clear that the amplitude 
of the stochastic fluctuations that act on the converted 
modes depend on the pump field dynamics. A brief dis¬ 
cussion of the noises is given in Appendix B. 


A. Critical driving field 


As a first investigation, we treat the classical approx- 
imation, which has also been analyzed in some earlier 


work [24 l3lL l54l - l57| . Here one assumes that all noise is 
negligible, so that At = A *, which gives equations in the 
form: 


dAo 

dt 

dA, 

dt 


~7oAq -I- £(x) — X*Ai A2 + —Aq , 

2a;o 

—liAi -l- xAoA3_j -|- —^V^Ai. (3-3) 

ZUJi 


The phases of £ and x are essentially arbitrary, as they 
depend on the phase definition for the field amplitudes 
Ai, which in turn depend on arbitrary mode phases. We 
will use this freedom later on to simplify the equations. 
If in addition, we assume that the input is a plane-wave, 
and we ignore possible spatial instabilities, diffraction can 
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be neglected as well. A steady-state result involves set¬ 
ting the time-derivatives to zero, so: 


A-o — (£ — x*AiA 2 ) /70 , 

A, = xAA*_Jl^ ■ (3.4) 


Hence, the defining equation for a steady-state is: 




^ 1^2 Ix^ol' 

7172 


(3.5) 


We can always choose the interaction picture detuning e 
so that 7 i 7 | = 7 ^ is real, otherwise the above threshold 
solutions will be oscillatory rather than stable. There 
are two types of steady-state solution. Either A 1 A 2 = 0, 
or else |xAo| = 7 . The first is called a below-threshold 
solution, the second an above-threshold solution. There 
is a driving field where both the solutions coincide, at a 
critical pump intensity of: 

l^cl^ = f l7o/xl^ • (3.6) 


More generally, suppose we are in the above-threshold 
regime. Using Eq. (EH), and defining = Ji, one 

obtains: 

Ao = (f-| x'|/iAo/72)/7o. (3.7) 

On re-arranging the equation, and using the above¬ 
threshold solution IxAqI = 7 , this result becomes: 


IxAqI^ =7^ = 


l72X'^^r 


I7072 + lx^| 7 i| 

The roots of the resulting quadratic are: 

1 


h = 


\X^ 




(3.8) 


(3.9) 


where z = z' + iz" = 7072 - 

Solutions with negative intensities are unphysical. 
There is a positive, above-threshold solution if \x£\ > \z\, 
which gives an identical critical field to Eq. (13.61) . Eor 
\£\ > £c there is a transfer of energy from the pump to 
the signal and idler modes, which develop a hnite mean 
intensity. It is the vicinity and just above this critical 
point that is the main regime of interest in this paper. 
We note that above threshold, further instabilities exist, 
including limit cycles and spatial pattern formation . 


IV. ADIABATIC ELIMINATION OF THE 
PUMP MODE 


is, we are neglecting pump diffraction), we can perform 
an adiabatic elimination by using the stationary solution 
for the pump mode, so that 


Aq = Aq = 


£--X*AiA2 

7o 


(4.1) 


A. Signal and idler equations 

The resulting equations for the down-converted modes 
- often called the signal and idler equations - are, for 
i = 1,2: 

A ' "v 

= - 7 ,A, + ^ (£■ - x*^i 7 l 2 ) At,+ 
ot 7o 

+xAQ^i{x,t). (4.2) 

We see that the main effect of detuning the pump is to 
reduce the effective intra-cavity pump intensity. So far, 
we have treated the case of general frequencies and group 
velocities. An important special case is obtained when 
the two down-converted frequencies are equal. In this 
case, the modes are still non-degenerate, as they can have 
different polarizations. From now on, we consider this 
special case for simplicity, so we take 71 = 72 = 7 , vi = 
V 2 = V and uji = UJ 2 = This implies that Ai = A 2 = 
A. We will also assume that Aq = 0, i.e. that the pump 
is on-resonance with the cavity, even when the down- 
converted helds may be off-resonant from their cavity 
resonance frequencies. 

Although these assumptions simplify the algebra, 
they are not essential for our main conclusions, which 
mostly rely simply on the fact that we now have a 
two-dimensional order parameter rather than a one¬ 
dimensional order parameter as found in the degenerate 
case. We should note that a non-zero pump detuning 
can excite another modes different from loq. This could 
give rise to nonlinear phenomena like bistabilities, non¬ 
linear resonances and subcritical bifurcation, especially 
for the case of large detuning . In the case where the 
diffraction terms are different, there will be a new term 
which is proportional to the difference of the two diffrac¬ 
tion terms, as has been studied elsewhere 0 . This term 
will be present in the equations even for the case of zero 
detuning. Since we are interesting in the universal be¬ 
haviour of the system we do not include these cases, but 
we point out that they can give rise to nonlinear be¬ 
haviour. 


We now return to the full quantum behavior given by 
the stochastic equations obtained above. One limit that 
has an especially simple behavior is found in the case of 
a rapidly decaying pump mode. We can treat this by 
means of an adiabatic elimination procedure. Assuming 
that 7 o ^ 7 i — 72 ) and that £ is spatially uniform (that 


B. Dimensionless form 

Equation (14.21) will now be transformed into a dimen¬ 
sionless form which allows comparisons with other types 
of phase transitions. First, we define the dimensionless 
variables r = t/tg, r = xjx^^ with a scaled Laplacian 
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/dr\ + ! dr 2 - It is useful to also define a 
dimensionless field = x^Ai as well. This has an in¬ 
tuitive interpretation as the coherent amplitude in real 
space, defined relative to a physical area of Xq. After this 
transformation, one obtains: 


1 dttj _ -r 

(4.3) 

We define 7 = 7 (1 + * A) and introduce the dimension¬ 
less pump amplitude 


A = — = . 

7o7 


(4.4) 


where /i is real and positive. Since the phase of the driv¬ 
ing field £ is arbitrary, and the equations are invariant 
under phase-changes of £, we will choose (j> = 0 with no 
loss of generality. We also define the time-scale to = 1/37 
as the scaling time for critical slowing down, where 


g^= 1^1' 

47072^0 

The length scale xq is now chosen so that: 


(4.5) 


27 V 5 W ■ 


(4.6) 


Combining Eqs. (I4.5|l and (14.6|) . one can write the di¬ 
mensionless coupling g in the form: 


9 = 



(4.7) 


Similarly, we introduce dimensionless complex noises = 
so that 

(Ci(’’,T)C 2 (r = 6'^{r - r ')S(t - r'), 

(Ci^ (r, t)C^ (r ', r')) = S^(r - r ')6{t - t'). (4.8) 

Finally, defining the driving field saturation factor as: 

/r (a) = ^ - 4g^aia2, 
fj,+ {d)= Ag'^at aj, (4.9) 


we obtain the following dimensionless form of the scaled 
equations for i = 1 , 2 : 


9^ = M («) “a'-i - (1 + 


+ V9 i^loLi A \/d {a)C,i{r,T) 


(4.10) 


together with a conjugate equation for . Although 
this is true in general, we are mostly interested here in 
the regime of g <C 1 , which allows us to make an expan¬ 
sion in powers of the coupling. For small g, the classical 
approximation gives the leading order term, while the 
diffraction and noise terms provide the next order in an 
expansion in 


V. STABILITY PROPERTIES AND 
QUADRATURE EQUATIONS 


We now wish to transform these equations into quadra¬ 
ture equations that are simpler to investigate. There 
are very different stability properties for the orthogonal 
quadratures near the critical point. To investigate this, 
as a first approximation, we will ignore noise and nonlin¬ 
ear terms of order and smaller. The stability of the 
equations near = 0 , to leading order in g, is: 



0.1 


(1 -I- zA) /i W “1 ^ 

9 - (1 - *A) j y ’ 

(5.1) 


which has eigenvalues A± = — 1 ± a/—A^ -|- /i^, and with 
a similar equation coupling 0.2 and a^. This gives an 
unstable eigenvalue, leading classically to growth of the 
signal and idler terms if — A^ > 1 , as expected from 
the analysis in the previous section. The resulting eigen¬ 
vectors, u±, are: 


u± = 


9 _\ 

zA± \/-A2 -h ) ■ 


(5.2) 


There is clearly a line of critical points where the stability 
has a continuous change at — A^ = 1. The tricritical 
point then occurs when /i = 1, A = 0, as we will show in 
later sections. 


A. Quadrature field variables 


To understand the behavior of Fq. (14.101) in the neigh¬ 
borhood of the critical point, we define complex, dimen¬ 
sionless scaled quadrature fields which are proportional 
to the critical eigenvectors, as follows l3^ : 


=\f9 (ai +at) , 

A+ =y/^{a2 +a^) , 

Y=i (ai -a+) , 

=- (a2 — a^) . (5.3) 

i 


Next, we consider the case ^ ft! 1 ^ 5 ^ 1 0102 !, to sim¬ 
plify the noise term. As this is both relatively small, and 
nearly constant in the neighborhood of the critical point, 
the resulting terms are of higher order in g than the lead¬ 
ing terms we wish to include. The resulting equations for 
these quadratures in the positive P-representation near 
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the critical point are: 

^ = —X + V+Y - (X^ + gY^) X+ + C+. 
or g 

= y x+ + V+Y+ - (x+2+ gY+^) X + c;. 

dY 

g— = -g+Y + V_X - g (X^ + gY^) Y+ - iy/gC-- 

9^ = - g (X+2 + gY+^) Y - z^C- 

(5.4) 

Here we have defined a modified Laplacian and driving 
term as: 


A , 
v± = ±—tvI 
V9 

g± = /i ± 1. (5.5) 

and new Gaussian noise terms according to: 

C±=Ci±C2+ = (C2±Ci+)*, (5.6) 

where we have used the result of Eq. (IC5I) . This shows 
that the noise terms driving the X, A+ fields are con¬ 
jugate, while those driving the Y, y+ fields change sign 
on conjugation, so that Y, Y^ will not remain conjugate 
during time-evolution. 

At this stage, we can make the following remarks. 
The stochastic quadrature fields X, A+ are both complex 
fields, so they have four degrees of freedom between them, 
and similarly for Y, T+. They have a correspondence 
with non-Hermitian operator fields X, and Y, Y^. 
In general, T, T+ are not complex conjugate except in 
the mean, and neither are X, A+, since they are driven 
by the Y, Y'^ fields. However, as we will show, this pic¬ 
ture simplifies when one considers an expansion near the 
critical point. 


B. Critical point adiabatic elimination 


We can now perform a second type of adiabatic elimi¬ 
nation which is valid in the neighborhood of the critical 
point. This takes into account the fact that the fluctua¬ 
tions in the X quadrature become very slow near thresh¬ 
old, while the Y quadrature still responds on the fast 
relative time scale I/ 7 . Formally, we can drop terms of 
0{yfg) where g <C 1, and approximate equations (15.41) as 
follows 


dX 

dr 

dX+ 

dr 



X + 
X+ + 


A 


V~9 



0 = - (1-H ^) r-f A, 

0 = -(l-h^) Y+ +Vlx+. 


Y - X^X+ + C+ , 

- Y+ - X+^X + c; , 


Here we have considered that the noise term of the 
quadrature variables Y, F+ has been neglected since it 
scales as y/g and we are considering the limit g ^ 1. If 
this limit is not taken, the noise term would appear in 
the above equations. To lowest order in g, we can elimi¬ 
nate the fast or non-critical quadrature Y, Y~^ variables 
by setting: 


y(+) 


1 1-1 


(5.8) 


which gives the result that: 

dX ~ , , 

— =VX- X^X+ + C+ , 

^=VX-{X+f X + C+- (5.9) 

Here we have introduced a linear differential operator 
that describes the linear gain, loss and diffraction terms: 

v = Vr = -vi + mXl-mxt (5.io) 


where we have defined the following parameters: 


m 

m 

m 



A 


{l+g)y/g' 

1 

1 + ^ 


(5.11) 


We note that since the non-conjugate variables T, F+ do 
not appear in these equations, it follows that variables 
X, X'^ will remain conjugate if they are conjugate ini¬ 
tially, as for example in an initial thermal or coherent 
state. Hence, to leading order, we can set A+ = X* near 
the critical point, and write one complex equation for the 
critical quadratures, which is valid near threshold: 

BX 

— =PA-A|A|2 + C+. (5.12) 

While this equation is similar to a Landau-Ginzburg 
equation for a complex order parameter [s^ . it is not 
identical, owing to the presence of the fourth-order Lapla¬ 
cian term va.'D. In the next section we will show that the 
equations are similar to a vector Swift-Hohenberg equa¬ 
tion [^. 


C. Vector Swift-Hohenberg equations 


We will now show that these near-threshold equations 
are actually coupled or vector Swift-Hohenberg equa¬ 
tions [ 2 | that represent the leading order dynamics near 
threshold of the down-converted modes with the same 
frequency but orthogonal polarization. To demonstrate 
this, it is convenient to make the following change of vari¬ 
ables: 


Ai 


A-f A* 


A 2 


A-A* 


(5.7) 


2 


2i 


(5.13) 














After this change of variables it is possible to write the 
equation (15.121) as a single vector equation, as we now 
show. On inverting this equation we get: 


X = Xi + iX2, X* =Xi- iX2. (5.14) 


Next we define define two real vectors X and C as: 


X = 




Using the above expressions we can write Eq. (15.1211 as 
a single vector equation of the form: 

= t>X-\X\^X + C (5.15) 

OT 

where X is a two real component vector whose elements 
are Xi and X 2 , and C is also a two real component Gaus¬ 
sian noise vector, with correlations given by: 


,t')^ = S^j6{r - r')S{T - t'). (5.16) 



Figure 1. Dimensionless intensity versus dimensionless pump 
in the vicinity of the critical point. The point p = 1 cor¬ 
responds to 771 = 772 = 0 and ris = ^. Here we have used 
the parameter g = 0.01. Results were obtained from a simu¬ 
lation with 300 samples on a 50 x 20 x 20 numerical grid of 
30000 X 96 X 96 points. 


In the case of a one-dimensional, real order parameter, 
this equation was first derived by Swift and Hohenberg [3, 
ilMl. Here it appears as a two-dimensional vector equa¬ 
tion, since the order parameter is two-dimensional. This 
was used to explain the convective roll patterns gener¬ 
ated by the Rayleigh-Benard instability, where the or¬ 
der parameter is the vertical fluid velocity. The real or¬ 
der parameter case is also similar to the Ising model for 
magnets in two dimensions with next nearest neighbor 
interactions where competition between the near¬ 
est and next nearest interactions generates a magnetic 
modulated phase called the Lifshitz phase (^ . 

Higher dimensional or complex order parameters, as 
in the present analysis, are described using a general¬ 
ized Landau-Ginzburg-Wilson Hamiltonian (^ . where 
the authors also introduce the Lifshitz point. For the 
present case, the upper critical dimension for classical 
behavior is at d = 5, and the classical location of the 
Lifshitz point is at 771 = 772 = 0. Low dimensional 
cases should have enhanced fluctuations, without sponta¬ 
neous magnetization or symmetry breaking, as expected 
from the Mermin-Wagner theorem (s^. This applies to 
Heisenberg-type models with higher-dimensional order 
parameters in two dimensions. As it is valid for general, 
finite-range interactions, it also holds in our case. 


VI. CORRELATIONS 

All physical quantities that we wish to understand in 
the double adiabatic limit treated in the previous section, 
come from the solution of equation (15.151) . This allows us 
to calculate the expectation value of any other observable 
in the vicinity of the critical point. Another way to obtain 
expectation values is write the functional probability as a 
solution of the master equation. Below we develop both 
methods. 


A. Stationary solution of the Fokker Planck 
equation 


In Fig. ([Ij we plot the modulus squared of the two- 
dimensional order parameter X vs the dimensionless 
pump parameter /r, averaged over a transverse area of 
dimensionless size 20 x 20. 

Direct simulations were employed, using both a cen¬ 
tral partial difference algorithm in the interaction pic¬ 
ture ( 5 ^ and a fourth order interaction picture Runge- 
Kutta method. Two public domain software packages 
were used and checked against each other [^. This fig¬ 
ure is obtained by solving the stochastic equations given 
in Eq. (15.151) at /r = 0.9, equilibrating for t = 50 units 
then scanning the driving field n adiabatically (using 
/j, = 0.9-|-0.004f) through the critical point until /r = 1.1. 
There is a rate-independent critical region at the thresh¬ 
old value of 77 = 1 , where the transition is smooth rather 
than discontinuous, as it would be classically. 


We start with equation (|5.15|) and note that it is possi¬ 
ble to write a functional Fokker-Planck equation for the 
probability density P{X,t), 


dP 


iX,r) 



[\x\^-v')x. 


1 6 
26 X 


P{X, 


( 6 . 1 ) 

and look for the equilibrium distribution in the form 
P{X) = Vexp[—'H(X)], where If is a potential func¬ 
tion. The solution for the distribution of X is given by: 


P{X) = N exp 


/ 


- Px( T]lX ■ X 


l{x-xf + 


+r]2VX ■ VX + rjsX^X ■ V^X)] . 


( 6 . 2 ) 


This expression is similar to the Landau-Ginzburg free 
energy of a next nearest neighbor interaction in a planar 
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magnetic interaction, with X playing the role of a two 
component vector order parameter. Owing to this par¬ 
allel, the planar type-II parametric system can provide a 
superb model platform for investigating fluctuations and 
universal behavior in this paradigmatic system. 


B. Stochastic moments in the Gaussian 
approximation 


In order to evaluate the moments and spatial corre¬ 
lations we will approximate the nonlinear terms of Eq. 
(|5.15l) using a Gaussian approximation together with 
a Green’s function approach. Hence we can write Eq. 
(|5.15l) as follows: 


d{X,{r,T)Xj{r',T)) 

dr 


T>r {Xi{r,T)Xj{r' ,t)) 

iXi{r,T) \X {r,T)\^ Xj{r',T)^ . 

(6.3) 


Here (...) denotes the stochastic average over the Gaus¬ 
sian fluctuations and we have used the notation 7)^ hi 
order to avoid ambiguities. We explicitly denote that 
the operator T) acts on the spatial coordinate r. In the 
Gaussian approximation, fluctuations around the most 
probable configuration are approximately treated as in¬ 
dependent modes with Gaussian distributions (^ . This 
allows us to replace ensemble averages by a Gaussian 
ansatz, in which higher order moments are approxi¬ 
mated by the expressions for a Gaussian distribution, 
e.g. {X'^) ~ 3{X'^)'^. On defining Gy = Gy {r,r') = 
{Xi{r,T)Xj{r',T)) we can write the above equation, for 
1 = 1 , as: 


= P,Gy - Gi, (3 (Xf) + (X|)) . 

Here we have assumed rotational symmetry of the prob¬ 
lem in the {Xi,X 2 ) plane, so that {X 1 X 2 ) = 0. Us¬ 
ing again the assumption of rotational symmetry in the 
Xi — X 2 plane, we can also write the Green’s function 
results as: 


dr 


f>r - 2 ((X^) + (X|)) 


G 


u- 


(6.4) 


This corresponds to an equivalent stochastic equation, 
valid in the steady-state for a rotationally symmetric sys¬ 
tem: 


= VX{r, r) - 2 (X • X) X(r, r) + C(r, r). 

(6.5) 

We use the variable X to denote that we are using the 
Gaussian approximation. In order to evaluate the Gaus¬ 
sian correlation functions in near and far field, we will 
define a parameter t][ = rji -|-2(X -X). Using the expres¬ 
sion for V of Eq. (15.1011 we can write the above equation 


as: 

( 6 . 6 ) 

Next, we note that because of translational symmetry the 
term ^X • X^ is independent of the spatial coordinate so 

that (^X ■X'j = 

In order to evaluate it, we use the Fourier transform 

X(r,t) = ^ J e^^-^Xik,t)d^k, (6.7) 

so that in momentum space we write Eq. dnn) as: 


X 


can be calculated analytically. 


=- (Vi+ mk'^ + rjsk^) X{k,T) + C{k,T). 

( 6 . 8 ) 


C. Lifshitz point 


We now consider the line of points where 772 = 0 and 
773 = i. These are the points we have defined as corre¬ 
sponding physically to zero detuning, with a pump in the 
vicinity of the fi = 1. In this case the solution for X(fc, t) 
is given by: 


X{k,T)= f C(fe,T')e + (6.9) 

J —00 


In this way, we obtain: 


(X.X) 


1 r°° 2 Trkdk 

47r2 Jo T][ + ^' 


( 6 . 10 ) 


On performing the integration there is a resulting self- 
consistency condition: 


(^X.x)= ^ ^ (6.II) 

402(771 + 2 (X • X)) 


At the Lifshitz point, where 771 = 0, we find that: 


X • X) = ( X ) = 0.25. 


( 6 . 12 ) 


As explained below, this is remarkably close to accurate 
numerical simulations of the correlations, including non- 
Gaussian fluctuations. 

We also consider the case where /r = 1 ,772 7 ^ 0 and 773 = 
i. This corresponds to the case of non-zero detuning 
A. In Fig. ([2) we plot the modulus squared of the two 
dimensional order parameter X vs the detuning A. 

In order to obtain this figure, we solve the stochas¬ 
tic equations given in Eq. (I5.15|) . and scan the detuning 
which is proportional to the parameter 772 defined in Eq. 
(15.111) . so that A = —0.5 -|- 0.005t. We notice that the 
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A 


Figure 2. Dimensionless intensity versus detuning A. The 
Lifshitz point corresponds to zero detuning. Here we have 
used /r = 1 and g = 0.01. Results were obtained from a 
simulation with 60 samples on a 200 x 50 x 50 numerical grid 
of 15000 X 96 X 96 points. 


fluctuations depends on the detuning. For a positive de¬ 
tuning there is a decrease of the fluctuations while for a 
negative detuning the fluctuations increase. A classical 
Swift-Hohenberg equation with a complex order parame¬ 
ter and nonzero detuning has been treated [gJ , as has the 
hyperbolic complex Swift-Hohenberg equationj^ and 
other related studies (^.1^. For negative detuning there 
is a ring with strong fluctuations, shown in figures ([3]) and 



Figure 3. Dimensionless intensity versus detuning A and 
transverse momentum k^. For negative detuning there is are 
fluctuations for lower momentum values. Results were ob¬ 
tained from a simulation with 800 samples on a 100 x 50 x 50 
numerical grid of 15000 x 96 x 96 points. Here we have used 
fj, — 1 and g = 0.01. 



Figure 4. Left: Dimensionless intensity versus transverse mo¬ 
mentum k for a fixed value of the detuning A = 0.5. Here 
we have used p = 1 and g = 0.01. A ring is formed for lower 
values of the momentum. The fluctuations (peaks) are not 
symmetric due to noise. Right: Dimensionless intensity ver¬ 
sus transverse momentum for a fixed value of the detuning 
A = 0.45, and ky = 0. Other parameters are as in Fig. 


D. Spatial correlations in the Gaussian 
approximation 


From Eq (16.91) . the Gaussian correlation function in 
the momentum space (far field), in the stationary regime 
is therefore given by: 


(X(k) X(k') ) 


1 J(k -b k') 

2 + ri^k* ' 


(6.13) 


This should be readily observable experimentally, since 
the far-field region of the output field from the planar cav¬ 
ity is simply the Fourier transform of the internal cavity 
field. The correlation function for the enhanced quadra¬ 
ture in configuration space, or near field, is: 


(X(r)X(r')) = ^ / d'k- 


o-*k-(r-r') 


Vi + ^2^^ + Vsk''^ 
On integrating over the angle variable we obtain: 


1 


(X(r)X(r')) = — 


dk 


kJo {k\r - r'\) 


47r Jo r][ -I- r]2k'^ -|- rj^k'^ ’ 


(6.14) 


(6.15) 


where Jq is the Bessel function of zero order. The result 
of the above integration is, in the limit 772 = 0 : 


(X(r)X(r')) = - 


V3 1 


Vi 47rr73 


kei 


-) 

mj 


1/4 


r — r 


_ (6.16) 
where kei(x) is Thomson’s function [ 6 g. An interesting 
remark is that as ri[ —>- 0 , the spatial correlation decays 
with a power law: 


(X(r)X(r')) oc |r — r' 


/|- 0.5 


(6.17) 


We note, however, that this system is predicted to have 
an upper critical dimension |28l | with mean-field critical 
exponents at spatial dimension d > 5. Since the sys¬ 
tem is well below this critical dimension, one may expect 
non-Gaussian behavior that is not predicted by the ap¬ 
proximations used in this section. 
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E. Non-Gaussian behavior and universality 


In order to verify these analytic results and investigate 
non-Gaussian correlations, numerical simulations were 
carried out of the original stochastic partial differential 
equations of Eq. (I5.15|l . Small time steps are needed to 
treat the quartic growth of the squared Laplacian term 
in momentum, together with large sample numbers to 
obtain a low sampling error. In initial investigations, we 
used alOxlOxlO numerical grid of 8000 x 64 x 64 
points in t, x, y respectively. Employing a fine numerical 
grid of 16000 x 64 x 64 points to check convergence in 
time-step, the final steady-state correlation result con¬ 
verged to {X ■ X) = 0.264 ± 0.005. This was close to 
the Gaussian value, but with a relatively large sampling 
error. 

In order to understand the quantitative difference be¬ 
tween the exact and Gaussian results, a more precise dif¬ 
ferencing technique was used. This variance reduction or 
differencing technique simulates the difference between 
the full sample path and the Gaussian approximation 
( 3 ^ . The results were in agreement with a direct simu¬ 
lation, but gave much more rapid convergence. This was 
carried out as follows. First a mean field variable, X was 
simulated, in the Gaussian approximation, using: 


dX 


= VX-2X 



- 2X* (X^) + C+ , 


(6.18) 


where the averages were carried out both over a spatial 
numerical grid and a set of parallel trajectories. From 
the analytic calculations above, this should converge pre¬ 
cisely to the analytic result of = 0.25 , which was 

confirmed to a numerical accuracy of ±5 x 10“^. Here 
we have considered that yi = r ]2 = 0 ^^nd rjs = ^, corre¬ 
sponding to the classical Lifshitz point. 

Next, the full variable X was simulated. This was 
achieved by introducing a difference variable defined as 
Ax = X — X, which has the stochastic equation: 



Figure 5. Growth of non-Gaussian correlations, {\Xf‘ — \X\'^) 
versus time r, starting from X = 0. Error bars were obtained 
from comparing a coarse (5000 step) and fine (10000 step) 
simulation. The sampling error is indicated by the upper and 
lower solid lines, with a standard deviation of ±0.00025. This 
gives the steady-state correlation result {X ■ X) = (|A|^) = 
0.2574 ± 0.0003. Results were obtained from a simulation 
with 3200 samples on a 10 x 20 x 20 fine numerical grid of 
10000 X 48 X 48 points. 




Figure 6. Left: Steady-state momentum correlations, 
(|A’(k)|^) versus transverse momentum k, starting from X — 
0. Right: Steady-state non-Gaussian correlations in momen¬ 
tum space, Ain (|A(k)|^) = In (|X(k)|^) — In ^|A’(k)|^^ ver¬ 
sus momentum k. Here we consider 771 = >72 = 0 and rjs — 
Other parameters are as in Fig. m 


^=VX-X\x\^ + C+-^. (6.19) 

By using an identical noise source to those in the equa¬ 
tions for the mean-field X, the difference simulation per¬ 
mits a more precise calculation with reduced variance. 
The result, shown in Figs. ([S]) and ([5]) was that, at the 
critical point, {X ■ X) = 0.2574 ± 0.0003. This result, of 
much greater accuracy, only required 3200 samples with 
a 10 X 20 X 20 numerical grid of 10000 x 48 x 48 points. 
The discretization error was estimated from using several 
grids with different transverse sizes and time-steps. 

In summary, the full statistical calculation gives in¬ 
creased critical fluctuations due to non-Gaussian effects, 
but this increase is relatively small. At large transverse 
momentum there is no measurable difference between the 
Gaussian and exact results, shown in Fig. ©• The devia¬ 
tion from the Gaussian approximation vanishes rapidly as 


higher order transverse momenta are investigated. Small 
values of this difference are observed only at zero trans¬ 
verse direction. 


VII. CONCLUSION 

We have shown that parametric down-conversion in a 
type-II parametric planar cavity leads to a Swift and Ho- 
henberg type of stochastic equation for the leading terms 
in the critical fluctuations, but with a vector order pa¬ 
rameter. This combines the rotationally invariant sym¬ 
metry properties of the X-Y model with the higher-order 
Laplacian of a Lifshitz magnetic phase transition. Sur¬ 
prisingly, these fluctuations are not thermal in origin, but 
come instead from the quantum fluctuations associated 
with parametric amplification. 
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This model can be approximately treated for the criti¬ 
cal fluctuations with a Gaussian factorization. However, 
a careful numerical treatment shows that non-Gaussian 
critical fluctuations occur. These are responsible for en¬ 
hanced intensity correlations, but are reduced at large 
transverse momenta due to the momentum dependence 
of the linear propagator. As techniques improve, we ex¬ 
pect that this novel, non-equilibrium critical point will 
become accessible to experimental studies. 

APPENDIX A: QUANTUM LANGEVIN FORM 

In the quantum Langevin form the corresponding op¬ 
erator equations of the system would be: 

= ~7o^o + £(x) — X* Ai A2 + V270^0" 

^ = -7iii + X^oAl + : 

^ = - 72^2 + X^4 + ^^'^2 + • 

(Al) 

Here we use a rotating frame such that the three field 
operators are treated as in a frame rotating with fre¬ 
quency uj^. The relative detuning between the pump 
laser at 2 ujl and the intracavity pumped mode ujq is 
Aq = (ojQ — 2 ujl) /701 and the down-converted modes 
have relative detunings Ai = (ui — /ji- The terms 
7 i = 7 i (1 + iAi) represent the complex cavity decay for 
each mode, including detunings. However, the quan¬ 
tum Langevin approach has the drawback that it deals 
with operator equations that are intractable analytically. 
It was shown in section m that the phase-space that 
the phase-space representation method generates similar 
equations, but with a more useful c-number form. 

The input-output relations that describe the external 
modes outside the cavity are d, [l^, = 

2'yiAi — A ■", where A™ and A°“* are the corresponding 
input and output fields, with input correlations: 

(^AY"{x,t)AY'\x\t)'j = (fil^ + l) 6ij6 (x - x '), 

= nfS,jS{x - x'). (A2) 

In our calculations we assume that the reservoirs are in 
the vacuum state. However, non-zero reservoir tempera¬ 
tures can be readily included. While we do not use these 
input-output equations here, we note that they are im¬ 
portant when dealing with external measurements. 

APPENDIX B: POSITIVE-P REPRESENTATION 

The positive P-representation generates a genuine (sec¬ 
ond order) Fokker-Planck equation with positive-definite 


diffusion, provided the distribution vanishes sufficiently 
rapidly at the phase-space boundaries. This can then 
be mapped into a set of c-number Langevin equations 
similar to the quantum Heisenberg equations, except for 
additional stochastic terms. 

This approach uses a multi-mode coherent state 
I Aq) 0 ^ 1 ; 0 : 2 ) = |q:), defined as an eigenstate of the an¬ 
nihilation operators ai (fe), where oti = oti (k) so that the 
following eigenvalue equation is obtained: 

at {k) \a) = Ui {k) \a) . (GI) 

The positive-P representation then expands the den¬ 
sity matrix in terms of coherent state projection oper¬ 
ators [13 , which is always possible as a positive distribu¬ 
tion: 

p = J ^ + ^ ^ 

In the positive P-representation the coherent ampli¬ 
tudes di (fc, t) satisfy time-dependent stochastic differ¬ 
ential equations. It is simplest to write these equations 
in a form analogous to classical equations by introducing 
stochastic fields Ai , A)*' defined as: 

A,(a:,t) = (G3) 

k 

together with a stochastic conjugate Af which is a c- 
number, rather than an operator field. It is only conju¬ 
gate to Ai in the mean: it is stochastically equivalent to 
the conjugate operator. 

1. Noises of the stochastic equations 

We note that while our derivation of the set of stochas¬ 
tic equations given in Eq. (Ell) is formally based 
on the Ito stochastic calculus, in this case either Ito 
or Stratonovich stochastic calculus gives identical re¬ 
sults. These complex noise terms can be constructed 
from four delta-correlated real Gaussian noise fields 
, Cjj"), with the mapping: 

6.2(a;A) = [(xix,t)±i^y{x,t)]/V2, 

= [i+{xA)±i^+{x,t)] IV 2 . (G4) 

It follows that the stochastic fields in the positive P- 
representation for the and ^2 fields are complex con¬ 
jugate, i.e., 

^i(a;,t) = , and {xA) = {it A))* ■ 

(G5) 

The physics of the noise is that it describes a departure 
from coherent behavior. The deterministic terms corre¬ 
spond to the evolution of coherent states. However, the 
true quantum state does not remain coherent. Instead it 
develops, via the noise terms, into squeezed, entangled, 
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and even more complex states. Despite this complexity, 
there are simple, universal properties caused by these 
noise terms at the critical point. 

We note that there are no ‘normal’ noise correlations, 
that is, {x',t')) = 0. This is due to our as¬ 

sumption that the optical reservoirs are at zero temper¬ 
ature, which is an excellent approximation at optical fre¬ 
quencies. If there are thermal reservoirs, as can occur in 
microwave devices, then additional reservoir correlations 
must be included, which are proportional to the thermal 
occupation number. More generally, our model includes 
only the minimal noise due to fundamental quantum ef¬ 


fects, but there can be a range of additional technical 
noise sources in practical devices, caused by temperature 
fluctuations, laser intensity fluctuations and laser phase 

noise 0113 . 
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